###Replication data for "Does inequality beget inequality? Experimental tests of the prediction that inequality increases system justification motivation"
### Kris-Stella Trump and Ariel White
### December 2017
### This replication code for Additional Study 1 presented in Supplemental Information
### Code tested in R version 3.4.1

#Create log file
sink(file="Add_Study1_analysis_code_log.txt")

#Packages
install.packages("psych")

require(psych)

#Loading data 
#Don't forget to set working directory
gdata = read.csv("SJT_Gini_data.csv")
names(gdata)

#### Basic descriptive stats:
# sample from this go-round:
gdata$yearborn <- 1995-gdata$born
summary(gdata$yearborn) 
table(gdata$gender) 
summary(gdata$educ) 
length(unique(gdata$state)) 
summary(gdata$hhsize) 



#Treatment variable recode
gdata$treat <- is.na(gdata$ginilow)
summary(gdata$treat)

#####Divide data by dependent variables 

####SJT outcome group
class(gdata$sjt_1)
summary(gdata$sjt_1)
sjt.gdata <- gdata[is.na(gdata$sjt_1)==F,]

#Reverse code items 3 and 7
sjt.gdata$sjt_3r <- 10-sjt.gdata$sjt_3
sjt.gdata$sjt_7r <- 10-sjt.gdata$sjt_7

sjt.gdata$sjt <- (sjt.gdata$sjt_1 + sjt.gdata$sjt_2 + sjt.gdata$sjt_3r + sjt.gdata$sjt_4 + sjt.gdata$sjt_5 + sjt.gdata$sjt_6 + sjt.gdata$sjt_7r +sjt.gdata$sjt_8)/8
summary(sjt.gdata$sjt)

sd(sjt.gdata$sjt, na.rm=T)
t.test(sjt.gdata[sjt.gdata$treat==0,"sjt"],sjt.gdata[sjt.gdata$treat==1,"sjt"])

#Internal consistency 
sjt <- sjt.gdata[c("sjt_1", "sjt_2", "sjt_3r", "sjt_4", "sjt_5", "sjt_6", "sjt_7r","sjt_8")]
alpha(sjt)

#####Gender outcome group
gender.gdata <- gdata[is.na(gdata$gender.inf)==F,]

#Gender dv - recode so higher numbers indicate more justification of differences
gender.gdata$gender.dv_1r <- 7-gender.gdata$gender.dv_1
gender.gdata$gender.dv_3r <- 7-gender.gdata$gender.dv_3
gender.gdata$gender.dv_4r <- 7-gender.gdata$gender.dv_4
gender.gdata$gender.dv_5r <- 7-gender.gdata$gender.dv_5

gender.gdata$gender <- (gender.gdata$gender.dv_1r + gender.gdata$gender.dv_2 + gender.gdata$gender.dv_3r + gender.gdata$gender.dv_4r + gender.gdata$gender.dv_5r)/5
summary(gender.gdata$gender)
sd(gender.gdata$gender, na.rm=T)
t.test(gender.gdata[gender.gdata$treat==0,"gender"],gender.gdata[gender.gdata$treat==1,"gender"])

#Internal consistency
gender <- gender.gdata[c("gender.dv_1r", "gender.dv_2", "gender.dv_3r", "gender.dv_4r","gender.dv_5r")]
alpha(gender)


######Trust outcome group
trust.gdata <- gdata[is.na(gdata$gss.conf_1)==F,]

#Trust DV - all in one direction, recode so higher number is more trust
trust.gdata$gss.conf <- (trust.gdata$gss.conf_1 + trust.gdata$gss.conf_2 + trust.gdata$gss.conf_3 + trust.gdata$gss.conf_4 + trust.gdata$gss.conf_5 + trust.gdata$gss.conf_6)/6
  trust.gdata$trust <- 4-trust.gdata$gss.conf

sd(trust.gdata$trust, na.rm=T)
t.test(trust.gdata[trust.gdata$treat==0,"trust"],trust.gdata[trust.gdata$treat==1,"trust"])

#Internal consistency
trust <- trust.gdata[c("gss.conf_1", "gss.conf_2", "gss.conf_3", "gss.conf_4", "gss.conf_5", "gss.conf_6")]
alpha(trust)

#Treatment groups
dim(sjt.gdata); dim(gender.gdata); dim(trust.gdata)
#85,86,85

###Results: putting t-tests from above all in one place:

t.test(sjt.gdata[sjt.gdata$treat==0,"sjt"],sjt.gdata[sjt.gdata$treat==1,"sjt"])
t.test(gender.gdata[gender.gdata$treat==0,"gender"],gender.gdata[gender.gdata$treat==1,"gender"])
t.test(trust.gdata[trust.gdata$treat==0,"trust"],trust.gdata[trust.gdata$treat==1,"trust"])

sd(sjt.gdata$sjt, na.rm=T)
mean(sjt.gdata$sjt, na.rm=T)

###Plotting results

#put everything on a comparable scale (compress to 0-1)
sjt.gdata$sjtnew <- sjt.gdata$sjt/9
gender.gdata$gendernew <- gender.gdata$gender/6
trust.gdata$trustnew <- trust.gdata$trust/4

##Look up the sd of each of the standardized measures. 
sd(sjt.gdata$sjtnew, na.rm=T)
sd(gender.gdata$gendernew, na.rm=T)
sd(trust.gdata$trustnew, na.rm=T)

##set up gini/various outcomes
gini_grid <- as.data.frame(cbind(rep(NA,3), rep(NA,3), rep(NA,3), rep(NA,3)))
gini_grid$V1 <- c("System \nJustification", "Gender \nJustification", "Institutional \nTrust")
colnames(gini_grid)<- c("outcome", "diff.est", "conf.low", "conf.high")

gini_sjt <- t.test(sjt.gdata[sjt.gdata$treat==1,"sjtnew"],sjt.gdata[sjt.gdata$treat==0,"sjtnew"])
gini_sjt$conf.int[1]
gini_grid[1, 2] <- gini_sjt$estimate[1]- gini_sjt$estimate[2]
gini_grid[1, 3] <- gini_sjt$conf.int[1]
gini_grid[1, 4] <- gini_sjt$conf.int[2]

gini_gen<- t.test(gender.gdata[gender.gdata$treat==1,"gendernew"],gender.gdata[gender.gdata$treat==0,"gendernew"])
gini_gen$conf.int[2]
gini_grid[2, 2] <- gini_gen$estimate[1]- gini_gen$estimate[2]
gini_grid[2, 3] <- gini_gen$conf.int[1]
gini_grid[2, 4] <- gini_gen$conf.int[2]

gini_trust <- t.test(trust.gdata[trust.gdata$treat==1,"trustnew"],trust.gdata[trust.gdata$treat==0,"trustnew"])
gini_grid[3, 2] <- gini_trust$estimate[1]- gini_trust$estimate[2]
gini_grid[3, 3] <- gini_trust$conf.int[1]
gini_grid[3, 4] <- gini_trust$conf.int[2]
gini_grid$ruler <- c(1,2,3)


##Figure A4 in Supplemental Information
pdf("gini_differenceests_std.pdf")
par(mar=c(5.1,6.1,4.1,2.1))
plot(NA,NA, main="Gini Treatment (Various outcomes)", ylim=c(.5, 3.5), xlim=c(-1,1), yaxt='n', 
	xlab= "Estimated Difference Between High/Low Treatment Conditions", ylab="")
points(gini_grid$diff.est, gini_grid$ruler, pch=20)
for (i in 1:3){
	segments(gini_grid$conf.low[i], gini_grid$ruler[i], gini_grid$conf.high[i], gini_grid$ruler[i])}
par(las=1)
axis(side=2, at=gini_grid$ruler, labels = gini_grid$outcome)
abline(v=0, lty=3, col="gray")
dev.off()


### time spent on treatments
names(gdata)
time <- c(gdata$ginilow.t_3, gdata$ginihigh.t_3)
summary(time) #most people spend definitely enough time
t.test(gdata$ginilow.t_3, gdata$ginihigh.t_3) #no difference in time spent by condition
sum(time>5, na.rm=T) #247 of 256 people spend more than 5 seconds


#### also, look at basic t-tests for just low-income people.

table(gdata$ind.inc1) #about 2/3 are below $35k/yr.
sum(gdata$ind.inc1==1)

#original t-tests/ new subsetted t-tests:
t.test(sjt.gdata[sjt.gdata$treat==0,"sjtnew"],sjt.gdata[sjt.gdata$treat==1,"sjtnew"])
t.test(sjt.gdata[sjt.gdata$treat==0 & sjt.gdata$ind.inc1 == 1,"sjtnew"],sjt.gdata[sjt.gdata$treat==1 & sjt.gdata$ind.inc1 == 1,"sjtnew"])

t.test(gender.gdata[gender.gdata$treat==0,"gendernew"],gender.gdata[gender.gdata$treat==1,"gendernew"])
t.test(gender.gdata[gender.gdata$treat==0 & gender.gdata$ind.inc1 == 1,"gendernew"],gender.gdata[gender.gdata$treat==1 & gender.gdata$ind.inc1 == 1,"gendernew"])

t.test(trust.gdata[trust.gdata$treat==0,"trustnew"],trust.gdata[trust.gdata$treat==1,"trustnew"])
t.test(trust.gdata[trust.gdata$treat==0 & trust.gdata$ind.inc1 == 1,"trustnew"],trust.gdata[trust.gdata$treat==1 & trust.gdata$ind.inc1 == 1,"trustnew"])
